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We present a mean field model for spin glasses with a natural notion of distance built in, namely, 
the Edwards- Anderson model on the diluted _D-dimensional unit hypercube in the limit of large 
D. We show that finite D effects are strongly dependent on the connectivity, being much smaller 
for a fixed coordination number. We solve the non trivial problem of generating these lattices. 
Afterwards, we numerically study the nonequilibrium dynamics of the mean field spin glass. Our 
three main findings are: (i) the dynamics is ruled by an infinite number of time-sectors, (ii) the 
aging dynamics consists on the growth of coherent domains with a non vanishing surface-volume 
ratio, and (iii) the propagator in Fourier space follows the law. We study as well finite D effects 
in the nonequilibrium dynamics, finding that a naive finite size scaling ansatz works surprisingly 
well. 



I. INTRODUCTION 

Spin Glasses (SG) are highly disordered magnetic sys- 
tems Rather than by their practical usefulness, SG 
are often studied as a paradigmatic example of a complex 
system. Indeed, they display an extremely slow dynamics 
on a complex free-energy landscape, with many degener- 
ate states. In addition, SG are a convenient experimental 
model for glassy behavior, due to the comparatively fast 
microscopic spin dynamics, as compared, for instance, 
with supercooled liquids. In fact, nowadays, SG are con- 
sidered as a playground to learn about general glassy 
behavior, minimization problems in Computing Science, 
biology, financial markets, etc. 

Maybe the most conspicuous feature of SG is Aging: 
they never reach thermal equilibrium in experimental 
times. Here we will only consider the simplest possible 
experimental protocol, the temperature quench (see 
for very interesting, more sophisticated experimental pro- 
cedures) : the sample is cooled below the critical temper- 
ature, Tc, and it is let to relax for a time at the work- 
ing temperature T. Its properties are studied at a later 
time t + t^. It turns out that, if a magnetic field was 
applied from the temperature quenched until tw, when 
it is switched off, the thermoremanent magnetization 
M{t,t^) decays with t/t^j (at least for 10"^ < t/t^ < 10 
and 50 s < < 10^ s 3]). This lack of a characteristic 
time beyond iw, the glassy system age, is known as Full 
Aging. We now know that Full Aging is an effective de- 
scription of the dynamics, no longer valid for t/t^ ^ 10'* 
and U ~ 10s 0. 

Nowadays, we know that the slow dynamics in SG is 
originated by a thermodynamic phase transition at Tc . 
Below Tc, the spins associate in coherent domains, whose 
size, ^(tw), grows with time. The lower the temperature, 
the slower the growth of ^(^w) is [experimentally, ^{t^ = 
100s;T = 0.9rc)~100 atomic spacings]. 

There is a lively theoretical debate on the proper- 
ties of the low temperature phase. Surprisingly enough, 
this controversy on the equilibrium properties of a non- 



accessible (in human time scales) spiiiglass phase is rel- 
evant to nonequilibrium experiments [6|. 

Mainly, there are two competing theories, the droplets 
0, and the replica symmetry breaking (RSB) one 
According to the droplets picture, the SG phase would 
be ferromagnetic-like, with a complicated spin texture, 
but essentially with only two equilibrium states. On the 
other hand, the RSB theory predicts an infinite number 
of degenerated states with an ultrametric organization. 
For both theories. Aging would be a coarsening process, 
in the sense that coherent domains of low temperature 
phase grow with time. The two theories disagree in their 
predictions for these domains properties. According to 
droplets, the domains would be compact objects, with a 
surface- volume ratio that vanishes in the high ^(^w) limit 
The SG order parameter is non zero inside of each 
domain. On the contrary, the RSB theory expects non- 
compact domains with a surface-volume ratio constant 
for large £,{t^). Furthermore, in a RSB system, the SG 
order parameter vanishes inside those domains. In re- 
cent times, a somehow intermediate theory (TNT), has 
been proposed [To], but, to our knowledge, no detailed 
dynamic predictions have been provided. 

The RSB theory is based on the Mean Field approxi- 
mation (MF) which, unlike in the ferromagnetic case, is 
highly non trivial. Indeed, after 30 years of study. Ag- 
ing is not yet fully quantitatively understood, not even 
within MF approximations. Therefore, non perturbative 
tools, such as Monte Carlo (MC) calculations, appear 
as an appealing alternative. Furthermore, MC calcula- 
tions are called for even at the MF level. Hence, one 
is interested in mean field models, that is to say models 
where the MF approximation becomes exact in the ther- 
modynamic limit (TL). The standard MF model, (the 
Sherrington-Kirkpatrick model, see [ll| and Sect.|TT|), has 
a number of disadvantages. It lacks a natural notion of 
distance (hence one cannot discuss a coherence length 
S,{tw)), or coordination number. Furthermore, its nu- 
merical simulation is computationally heavier. In fact, 
recent advances on the analytical study of spin glasses 
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on Bethe lattices has shifted the attention to these 
far more numericaUy tractable models which share with 
experimental systems the notion of a coordination num- 
ber. 

Here we wish to present a new MF model for spin- 
glasses: the spin-glass on a D-dimensional hypercube 
with fixed connectivity [l3l | . In the thermodynamic limit 
(which coincides with the large D limit for this model), 
Bethe approximation becomes exact. As a consequence, 
the statics is of Bethe-lattice type and can be computed. 
A nice feature of this new model, is that it has a natural 
definition of distance, which allow us to study spatial cor- 
relations within MF approximation. In other words: this 
MF model is more similar to a real D = 3 system than any 
other studied before or, at least, than those considered 
previously, since the space-time correlation functions can 
be studied. 

The structure of this paper will be the following. In 
Sect. |lT]we will describe the model and compare it with 
other MF models. In particular, in Sect. Ill Al we will 
address the problem of fixing the connectivity in a di- 
luted hypercube. In Sect. IIIII we will briefly explain the 
numerical methods we have used and in Sect. llVl we will 
introduce the observables measured during the simula- 
tions. In Sect. rVl and IVll the numerical results will be 
presented, both in equilibrium (as a test of the model) 
and nonequilibrium, respectively. In Sect. IVIII we will 
discuss finite size effects. The analysis will reveal the 
propagator jl4i |. for the first time in a numerical investi- 
gation. Our conclusions will be presented in Sect. I Villi 
Finally, we include extended discussions in two appen- 
dices. 



II. MODELS 

The standard model in SG is the Edwards-Anderson 
(EA) model. We will work with two kinds of degrees of 
freedom: dynamical and quenched. The dynamical ones 
correspond to the spins, ai, with z = 1, 2, . . . A^. We will 
consider them as Ising variables ±1. The non dynamical 
(or quenched) represent the material impurities. We will 
consider here two types of them: the connectivity matrix, 
nik — riki = 1,0 {fiik = 1 as long as spins i and k interact), 
and the coupling constants, Jik = Jki, which shall take 
only two opposite values (in general with certain excep- 
tions that will be discussed in the next paragraph, we 
will consider Jik = ±1, which defines our energy scale). 
The interaction energy is 

H = -'^^Jiknik'Jiak ■ (1) 

i<k 

Since the impurity diffusion time is huge compared to 
the spin- flip (picosecond), we will always work within 
the so-called quenched approximation: spins cannot have 
any kind of influence over the material impurities. Then, 
both the set of coupling constants in the Hamiltonian 
and its associated Gibbs free energy, will be considered 



random variables. Therefore, in order to rationalize the 
experiments, the useful free energy will be an average 
over the disorder. We will refer to each assignment of 
{Tiiki Jik} as a sample. Its probability distribution defines 
the actual EA model. The average over samples will be 
represented as (...). 

Only a few exact results are known for the Hamiltonian 
(H)), and all them were obtained within the MF approxi- 
mation Q . This approximation becomes exact in a weak 
infinite-range interaction model (in a ferromagnet n^k = 1 
and Jik = 1/A^ for each couple i, k). In SG it is usually 
represented as Uik = 1 for every couple i, k and, because 
of the random ferromagnetic and antiferromagnetic inter- 
action character, {Jik} are independent gaussian random 
variables, with zero mean and variance 1/N 

Computer simulations of long-range models are ex- 
tremely hard, because the energy evaluation for a system 
of N spins requires A''^ operations. The situation has 
improved since the discovery that EA models on Bethe 
lattices (not to be confused with Bethe trees) undergo 
Replica Symmetry Breaking at [T^ . 

A very popular realization of a Bethe-lattice spin glass 
is the EA model on a Poisson graph. A simple way 
of drawing one of these graphs consists on connecting 
{riik = 1) each possible pair of spins, i,k, (there are 
N{N — l)/2 possible couples), with probability z/(N~l). 
Thus, the number of neighbors of spin i, its coordination 
number follows in the large- A'^ limit a Poisson distribu- 
tion function with average z (the connectivity). We will 
consider z = 6 to mimic a three dimensional system. This 
kind of graphs are locally cycle- less: the mean shortest 
length among all the closed loops that passes through a 
given point is 0(logA''), i.e. the system is still locally 
tree-like. Computationally convenient as they are, Pois- 
son graphs still lack a natural notion of distance. 

A simple alternative consists on formulating the model 
on a Z3-dimensional unit hypercube. Thus, the spins are 
located in each of the hypercube vertices (then, N = 2^) 
and the bonds lie on the edges. Therefore, each spin can 
be connected with, at most, D = \og2 N spins. By anal- 
ogy with the Poissonian graph, we consider that a link is 
active (i.e. Uik — 1) over each edge with probability z/D. 
We call this model random connectivity hypercube. It is 
easy to prove that it is locally tree-like as well: the den- 
sity of closed loops of length I decays, at least, with D~'^ 
(i.e. with the squared logarithm of N, as it also happens 
in the Poisson graph). Incidentally, one could consider 
as well a non diluted hypercube, but this would have two 
shortcomings: the connection with three dimensional sys- 
tems would get lost, and the numerical simulation would 
become computationally heavy for large D. 

Note that, at variance with other infinite-dimensional 
graphs, the hypercube has at least two natural no- 
tions of distance: Euclidean metrics and the postman 
metrics [Slj. The two distances are essentially equivalent, 
since the Euclidean distance between two sites in the hy- 
percube is merely the square root of the postman dis- 
tance. In the following we shall use the postman metrics. 



which has some amusing consequences. For instance, our 
correlation- length will be the square of the Euclidean one, 
thus yielding a critical exponent v — I, doubling the ex- 
pected i^MF = 1/2. Of course, if we use the Euclidean 
metric we recover the usual exponent — 1/2. 

However, it turns out that the random connectivity 
hypercube suffers a major disadvantage. The inverse of 
the critical temperature in a ferromagnet [15| or in a SG 
[l6l | can be computed within the Bethe approximation: 



Kf"^ = atanh 
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, K^'^ = atanh. 
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(2) 



In this expression (n) ^ is a conditional expectation value 
for n, the coordination number of a given site in the 
graph. This conditional expectation value is computed 
knowing for sure that our site is connected to another 
specific site (this is different from the average number of 
neighbors of a site that has at least one neighbor!). A 
simple calculation shows that (n)-^ = 1 + z — in the 
random connectivity model. Since D =^ log2 N, we must 
expect huge finite size corrections (0(1/ log A^)) at the 
critical point. Note that this problem is far less dramatic 
for a Poisson graph where (n) ^ — 1 + z — . 

The cure seems rather obvious: place the occupied 
links in the hypercube in such a way that n = z (here, 
z = 6). Unfortunately, drawing these graphs poses a non 
trivial problem in Computer Science [17|. Our solution 
to this problem is discussed in the next paragraph. 



A. The fixed connectivity hypercube 

We have not found any systematic way of activating 
links in the hypercube that respects the fixed connec- 
tivity condition. Thus, we have adopted an operational 
approach: the distribution of bonds is obtained by means 
of a dynamic MC. We must define a MC procedure that 
generates a set of graphs that remains invariant under all 
symmetry transformations of the hypercube group. 

Specifically, we start with an initial condition in which 
all bonds along the directions 1 to 6 are activated (of 
course, this procedure makes sense only for D > 6). 
Clearly enough, the initial condition verifies the con- 
straint n = 6. We shall modify the bond distribution by 
means of movements that do not change n. We perform 
what we called a "plaquette" transformation (a plaquette 
is the shortest possible loop in the hypercube, of length 
4). We randomly pick, with uniform probability, one hy- 
percube plaquette. In case this plaquette contains only 
two parallel active links {uik = 1), these two links are 
deactivated at the same time that the other two are ac- 
tivated. On the opposite case, nothing is done [s^]. This 
guarantees that the set of generated graphs is isotropic. 

In order to this procedure to be useful, the dynamic 
MC correlations times must be short. In Fig. [H we show 
the MC evolution of the system isotropy. We make fcA^ 
plaquette transformations, and we control the density of 
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Figure 1: Generation algorithm of fixed connectivity graphs: 
for two system sizes [D = 20, 22) and two spatial directions 
(u = 1, 7), we represent the density of occupied edges as func- 
tion of the MC time. As MC time goes on, the system recover 
the lost isotropy induced by the initial condition. 



occupied bonds in two directions: the first direction (ini- 
tially occupied in every vertex) and the seventh direction 
(initially unoccupied). As we see, for two different sys- 
tem sizes, we get short isotropization exponential times 
(for I? = 22 we get Texp ~ 4.7A). For this reason, we 
assume that taking A: = 100 is long enough to ensure that 
the configurations obtained are completely independent 
from the initial condition. 

At this point, the question arises of the completeness 
of the generated set of graphs. We first note that our set 
contains proper subsets that are also isotropic. How- 
ever, finding them would require more involved algo- 
rithms which will not pay in a reduction of statistical 
errors (as we will show below, most of the sample dis- 
persion is induced by the coupling matrix {Jik\)- On 
the other hand, one could think that there are lacking 
graphs in our algorithm for a simple reason. The plaque- 
tte transformation cannot break loops: when we inter- 
change neighboring links we can only either join two dif- 
ferent loops or split up a loop into two loops. Due to the 
hypercube boundary conditions, in the initial configura- 
tion all sites belonged to closed loops. This situation can- 
not be changed by plaquette transformations. However, 
this objection does not resist a close inspection. In fact, 
a non-closed lattice path formed by occupied links should 
have an ending point with an odd coordination number, 
which violates the constraint n = z for any even z. Thus, 
all lattice paths compatible with our fixed connectivity 
constraint, do form closed loops. This argument, as well 
as the numerical checks reported below, make us confi- 
dent that the set of generated graphs is general enough 
for our purposes. Actually, we conjecture that our al- 
gorithm generates all possible fixed connectivity graphs 
with z even. 

One may worry as well about the applicability of the 
Bethe approximation to the fixed connectivity model, 
since all loops are closed. Actually, the crucial point to 
apply the Bethe approximation is that the probability of 
having a closed path of any fixed length should vanish in 
the large D limit. This is easy to prove for the random 
connectivity model. In the fixed connectivity case, one 
may argue as follows. Let us imagine a walk over the 
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Figure 2: For (top) random-z graphs and (bottom) Fixed- 
2 for plaquette transformations, the probability distribution 
function of the length of the second shortest path joining 
nearest-neighbors for different sets of graphs (mind the ver- 
tical axis is in logarithmic scale). Lines has been slightly 
displaced in order to help the visualization. 
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Figure 3: Comparison of finite volume effects in the crit- 
ical point estimators for ferromagnetic Ising model, 
both in the random (red dots) and fixed (green 
crosses)connectivity hypercubes. As a guide to the eye, we 
have included two different scalings with D. The estimator 
corresponds to the average of the inverse temperatures 
at which the Binder cumulant, Eq. ([7]), reaches the values 1.2 
and 2.4. 



nectivity hypercubes instead, where the finite volume ef- 
fects are only caused by the residual presence of short 
closed loops. 



closed path. On the very first step, the probability that 
the chosen link is present is z/D, whereas in the following 
step the probability of finding the link is (z — 1)/(_D — 1) 
in the limit of large D (since one of the z links available 
at the present site was already used to get there). This 
estimate implicitly assumes that the occupancy of differ- 
ent links is statistically independent. The independency 
approximately holds for large D and becomes exact in 
the D —> oo limit, where occupied links form a diluted 
set. At this point, the estimate of the number of paths 
of any given fixed length in the large D limit can be per- 
formed as in the random-connectivity case. One finds as 
well, in the fixed connectivity case, that the number of 
closed loops per site of a given length decays at least as 
0{1/D^). 

In addition to the above considerations, we may nu- 
merically compute in our graphs the length of the second 
shortest path that joins two connected nearest neighbors 
in the hypercube. In Fig. [2l we compare the probabili- 
ties for the length of such paths in the random (top) and 
fixed (bottom) connectivity models, for different system 
sizes [33|. In both cases, we note that the maximum 
of the probability shifts to larger length as D grows. We 
note as well that, for fixed connectivity, no tree- like graph 
arises [3^. 

A summary of our efforts is shown in Fig. [31 where 
we plot the evolution of the critical point with D for the 
ferromagnetic Ising model, defined on hypercubes with 
both random and fixed connectivity. As expected, the 
random connectivity model suffers very important finite 
volume corrections which make it essentially useless for 
numerical studies. The problem is solved using fixed con- 



III. NUMERICAL METHODS 

We have simulated the Hamiltonian ^ using a 
Metropolis algorithm ^18i] . In addition, we use Multispin 
Coding: since spins are binary variables, we can simulta- 
neously codify 64 systems in one single 64 bits word (all 
of them share the same connectivity matrix, Uik). With 
this common matrix, we find errors which are ~7 times 
smaller than those obtained with one single sample per 
matrix. This should be compared with the factor 8 that 
we would obtain in case no correlation were induced. Our 
program needs 0.29 ns/ spin- flip in an Intel i7 at 2.93GHz 
(in Ref. they report '^1.2 ns/ spin- flip on an Opteron 
at 2.0 GHz, for the simulation of the D ~ 3 EA model in 
the cubic lattice). 

In a nonequilibrium dynamical study such as ours, one 
computes both one-time and two-times quantities, see 
Sect. IIVI The calculation of two-times quantities implies 
the storage on disk of intermediate configurations. Disk 
capacity turned out to be the main limiting factor for 
the simulation. For this reason, we have worked in par- 
allel with two program versions: one valid for measuring 
quantities at one and two times and another restricted 
to the computation of one-time quantities. 

We have computed two-time quantities at temperature 
T = 0.7Tc, on systems with D = 16, 18, 20 and 22. The 
number of simulated samples were 8x64 samples for each 
system size (hence, for self-averaging quantities the sta- 
tistical quality of our data grow with D). 

Besides, since this new model requires intensive test- 
ing, we have computed equilibrium one-time quantities at 
T/Tc = 0.95,0.97,0.99,1,1.1,1.2,1.3 and 1.4. The sys- 



5 



tem sizes were again D = 16, 18, 20 and 22. The number 
of simulated samples was 128x64 samples per tempera- 
ture (at Tc we computed 256x64 samples). 



IV. OBSERVABLES 

The Hamiltonian ([T]) has a global symmetry Z2 (cTi 
—Gi for all i). Not as obvious is the gauge symmetry 
induced by the average over couplings. In fact, choosing 
randomly a sign for each position, Si = ±1, the energy 
(fT|) is invariant under the transformation 



J, 



ik 



ik • 



(3) 



Now, since the transformed couphngs SiSkJik are just as 
probable as the original ones, we need to define observ- 
ables that are invariant under the gauge transformation 
With this aim we form gauge invariant fields from 
two systems at equal time, that evolve independently 
with the same couplings, {cr|^\o-|^^} (real rephcas) or, 
alternatively, from a single system considered at two dif- 
ferent times: 



(4) 



We can define three kinds of quantities with both fields. 
1. One-time-quantities. The order parameter 



(5) 



vanishes in the nonequilibrium regime (the system is 
much bigger than the coherence length, ^(iw))- The non 
linear susceptibility is proportional to the SG suscepti- 
bility: 



(6) 



that grows with a power of £,{t^). The Binder parameter 
provide us with information about the fluctuations 



(7) 



In the Gaussian regime B — 3. In a ferromagnetic phase, 
B = I. In the SG phase, in equilibrium (that for finite 
volume corresponds to t^, — * 00), B grows with the tem- 
perature from B = 1 at T = 0. The equilibrium param- 
agnetic phase is in Gaussian regime. 

2. Two-time-quantities. The correlation spin func- 
tion tells us about the memory kept by the system, at 
time t + tvj, about the configuration at iw [IE] '■ 



(8) 
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Figure 4: C4(r, tw), Eq. (|lip . for tv, — 2^ and different system 
sizes, N = 2°, at T = Q.7Tc. 

The link correlation function 

Clink(i, <w) ^w)Cfc(t, iw) , (9) 

ik 

carries the information of the density of the interfaces 
between coherent domains at t^, that at t -\- t^, have 
flipped. In case the surface-volume ratio decayed with 
a negative power of £(tw) (droplets), Cii„k would be- 
come i- independent [2^. On the contrary, in a RSB 
system, Ciink = a + bC'^. Note that, for the Sherrington- 
Kirkpatrick model, one trivially flnds Ciink = C^, but the 
linear relation is not straight-forward in flxed connectiv- 
ity mean field models. 

3. Spatial correlation functions. In the unit 
hypercube, the binary decomposition of the spin index 
i — 1,2, ... ,2^ can be equal to its Euclidean coordinates, 
X. The spatial correlation function is 



04(7-, tw) = -^y^gx(^w)gx+r(<w) ■ 



(10) 



We consider r = \r\ as the distance in the postman 
metrics. It would look rather natural to average all 
the Ci{r,tyj) over the Nr — 
r=\r\: 



( J. ) displacements of length 



C4(r,tw) = -^ C4(r,iw). 



(11) 



However, see Fig. m C4{r, t^) does not present a hmiting 
behavior with D for a given t^. 

We can get a clue by looking at xsg(^w), Fig- El 
which does reach a thermodynamic limit. Since Xsg(^w) 
is nothing but the integral of Ci{r, t^) with the Jacobian 
(^) , it seems reasonable to define the spatial correlation 
function instead: 



E 

r. \r\—r 



Ci{r,U 



(12) 



The susceptibility is x(tj = 27r/t,tw) oc [1 - C{t,t^)]/T. 
On the other hand, when is fixed, C{t,t^) is just the 
thermoremanent magnetization [36|. 



We can see that Ci{r, t^) does reach the high-D limit. 
Figure [H at least for short ty,. Besides, in the para- 
magnetic phase, it is possible to compute analytically 
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Figure 5: SG susceptibility a T = 0.7Tc as function of t„ for 
different system sizes, N = 2^. 
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2.44(3) 
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Table 1: Comparison between the SG susceptibility in large 
D limit for the paramagnetic phase, Eq. ([14]), and numerical 
results for D = 20, 22. 



V. EQUILIBRIUM RESULTS 

Since the present work is the first study ever made 
of a EA model on a fixed connectivity hypercube it is 
necessary to make a few consistency checks. Equihbrium 
results are most convenient in this respect, since we have 
analytical computations (valid only for the large D limit) 
to compare with. 

We will briefly study the spatial correlations in the 
paramagnetic phase. In addition, we will check, by ap- 
proaching to Tc from the SG phase, that the SG transi- 
tion does lie on the predicted Tc, Eq. 



Figure 6: C4(r, t„), Eq. Jill), for D = 10 and 22 for = 2*, 2** 
and 2^^ at T = 0.7rc. This has to be compared with the 
behavior of C4(r', tw), Fig. |4] 



C4(r, see Appendix [XI taking first the limit ^ oo 
and making afterwards D oo. The resulting corre- 
lation function, which is only valid in the paramagnetic 
phase, is a simple exponential. Hence, both the equilib- 
rium and the nonequilibrium computations, suggest that 
one should focus on C4 rather than on C4. 

We note in Fig. [51 that in the SG phase, C4 is non 
monotonically decreasing with r, but rather presents a 
maximum. This maximum moves to bigger r with t^, 
then, the system has a characteristic length that increases 
with time. Thus, we can estimate the coherence length, 
by means of the integral estimator ^0,1 (^w): 



Co,i(iw) 



/g°°dr r C4(r,tw) 

W J 



(13) 



A major advantage of ^0,1 over more heuristic definitions 
of the coherence length, is that it is computed from self- 
averaging quantities (see details in [IS, HH, we note that, 
in this work, we have not tried to estimate the contribu- 
tion to the integrals by the noise-induced long distance 
cutoff) . 

The existence of such a characteristic length is the 
main advantage of the hypercube model over other MF 
models. 



A. Paramagnetic Phase 

Our very first check will be the comparison between 
the Monte Carlo estimate of the SG susceptibility (at 
finite D) with the analytical computation for infinite D: 



X{T) = 1 + 



; tanh T 



2 rp-l 



(14) 



1 - (z- l)tanh^T-i' 

see Appendix [A] Our results are presented in Table Hi 
We see that finite size effects increase while approaching 
Tc- For our larger system, D = 22, the susceptibility 
significantly deviates from the asymptotic result only in 
the range Tc<T < l.lTc. 

After the fast convergence to the large D limit observed 
in the SG susceptibility, the results for C4 are a little bit 
disappointing. In Fig. [Tjwe display C4(r, D)—Ci{r, 00) as 
a function of r. We can see that finite size effects become 
more important once one approaches Tc. Besides, finite 
D corrections as a function of r oscillate between positive 
and negative values. This is not surprising: the finite D 
corrections to the susceptibility, which are very small, are 
just the integral under these curves. More quantitatively, 
we see in Table HIl that the corrections with D for r — 1,2 
are 0(0"^). Indeed, the path counting arguments in 
Appendix [K\ are plagued by corrections of 0{D^^). 



B. SG phase 

In the SG phase, our test has been restricted to a check 
of Eq. ([2), that predicts a SG phase transition for the 
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Figure 7: Difference between the numerical and analytical 
spacial correlation function for different system sizes at two 
temperatures T = l.lTc and T = lAT^. 



VI. NONEQUILIBRIUM RESULTS 

In this section we will address the main features of the 
nonequihbrium dynamics obtained in our largest system, 
D — 22. The issue of finite D corrections will be post- 
poned to Sect, rvm 



A. The structure of isothermal aging 

The picture of isothermal aging dynamics in MF mod- 
els of SG behavior was largely drawn in (see also 
[23j). The dynamics is ruled by an infinite number of 
time- sectors: 



(15) 





r 


= 1 


r 


= 2 


D 


T = l.lTc 


T = 1.4Tc 


T = l.lTc 


T = 1.4rc 


16 


0.783(6) 


0.198(5) 


2.130(18) 


0.320(12) 


18 


0.779(4) 


0.201(3) 


2.115(11) 


0.327(7) 


20 


0.784(2) 


0.202(2) 


2.109(6) 


0.332(4) 


22 


0.7776(12) 


0.2006(9) 


2.083(4) 


0.324(2) 



Table II: D times the difference between C4(r), for finite D 
and infinite D, as computed for r — 1,2. The absence of any 
D evolution evidences finite-D corrections of order 1/D. 



high-D limit. With this aim, we compute the Binder 
cumulant, B{T), nearby Tc. For all T < Tc, we expect 
B{T) <3 for large enough D. As we show in Fig. [51 B{T) 
decreases with T and shows sizeable finite size effects. 
In fact, at r = 0.99Tc, we need to simulate lattices as 
large as 13 = 20 to find values below 3. Right at Tc, the 
Gaussian value B{T) = 3 is found for all the simulated 
sizes. 
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Figure 8: Equilibrium values of the Binder cumulant, Eq. 
(0, for several system sizes, as function of the temperature 
in units of the exact asymptotic value of Tc, Eq. 0, in the 
SG phase. 



The scaling functions /; are positive, monotonically de- 
creasing and normalized, i.e. 1 = The un- 
specified functions hi are such that, in the large limit, 
hi{t^)/hi{t + tyj) is 1 if t ^ i^' , while it tends to zero if 
t^tl^^. In other words, the decay of C between values 
Ci and Ci+i is ruled by the scaling function fi and takes 
place in the time-sector t ^ . 

This picture is radically different to the Full Aging of- 
ten found both in experiments and in 3D simulations. A 
full aging dynamics is ruled only by two sectors of time, 
Hi = 0,1. Nevertheless, recent experimental studies Q 
show that full aging is no longer fulfilled for t^ty^. Prob- 
ably more time-sectors must be considered to rationalize 
these experiments. 

However, Eq. (|15p is probably an oversimplification, 
since the spectrum of exponents /ii might be continuous. 
An explicit realization of this idea was found in the criti- 
cal dynamics of the trap model [l^l , where the correlation 
function behaves for large t^ as 

C(i,tw) = /(a(i,tw)) , a(t,iw) = logt/logtw. (16) 

Again, the scaling function / is positive and monotoni- 
cally decreasing. Clearly enough, in the limit of large t^ 
and for any positive exponent ^, if i = Ai(^, the corre- 
lation function takes a value that depends only on fi, no 
matter the value of the amplitude A. 

As expected, C{t,t^) is clearly not a function of t/t^ 
in our model, see Fig. [^1 On the contrary, data seem 
to tend to a constant value when — > oo in any finite 
range of the variable t/t^. This is precisely what one 
would expect in a time-sectors scheme. On the other 
hand, if we try (without any supporting argument) the 
Bcrtin-Bouchaud scaling, Eq. p6)) . see Fig. [TOl the data 
collapse is surprisingly good. Therefore, the nonequihb- 
rium dynamics in the SG phase seems ruled by a, not 
only infinite but continuous, spectrum of time-sectors. 

We note en passant that the scaling (|16p is ultrametric 
only if the scaling function vanishes for all a{t,t^) > 1, 
for details see Appendix IB] In fact, dynamic ultrametric- 
ity is a geometric property [l^l that states that for all 
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Figure 9: C{t,t^) over t/t^ for D = 22 and T = Q.7Tc. 
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Figure 10: Same data of Fig. [O] as a function of ait^t^), 
defined in Eq. (|16p . The window is a zoomed image of the 
central region. 



triplet of times t\ ^ ^2 3> ts, one has in the hniit ts — > cxd: 

^?(tl-t3,t3) =niin{C(ti -t2,<2),C(t2 -^3,^3)}- (17) 

Finding dynamical ultrametricity in concrete models has 
been rather elusive up to now. An outstanding example 
is the critical trap model (l^l, where /(a > 1) = 0. It 
is amusing that the trap model is not ultrametric from 
the point of view of the equilibrium states [l^. Thus, 
the casual connections between static and dynamic ul- 
trametricity are unclear to us. At any rate, since our 
scaling function in Fig. [10] does not show any tendency 
to vanish for oiit^t^) > 1, we do not find compelling 
evidences for dynamic ultrametricity in this model. 

We have also looked directly to the plots of C(ti— 12, ^2) 
versus C{t2 — ^3,^3) (see Appendix |B|) and we have not 
found convincing indications for the onset of dynami- 
cal ultrametricity. In this respect, it is worth to re- 
call similarly inconclusive numerical investigations of the 
Sherrington-Kirkpatrick model [26| . There are two pos- 
sible conclusions: 

1. the model does not satisfy dynamical ultrametricity 
in spite of the fact that it satisfies (according to the 
standard wisdom) static ultrametricity. 
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Figure 11: Cunk over C^(t,tw) for different tw at T = O.TTc 
and for _D = 22. 



2. Dynamical ultrametricity holds but its onset is ter- 
rible slow. 

Both conclusions imply that it is rather difficult to use 
the dynamic experimental data (or any kind of data) to 
get conclusions on static ultrametricy. Of course it would 
be crucial to check if static ultrametricity is satisfied in 
this model, but this is beyond the scope of this paper. 



B. Aging in Cunk 

Just as in the 3D case [2l|, the aging dynamics in SG 
in the hypercube is a domain-growth process, see Fig llTI 
For any such process, the question of the ratio surface- 
volume arises. When this ratio vanishes in the limit of 
large domain size, as it is the case for any RSB dynam- 
ics, one expects a linear relation between Ciink(t, tw) and 
C^(t, iw). This is precisely what we find in Fig. [TT] 



C. Thermoremanent magnetization 

The experimental work indicates that for T < 0.9Tc, 
the thermoremanent magnetization follows a power law 
with an exponent proportional to T^jT [l^. The data 
obtained in JANUS for a three dimensional SG (see Fig. 
[T2land [20^) agree with this statement. However, the data 
obtained in the hypercube model does not follow such 
power law, neither can them be rescaled with T log t. 

This lack of an algebraic decay is surprising on the view 
of the exact results of Ref. [2^. Indeed, it was analyti- 
cally shown there that, at Tc, the thermoremanent mag- 
netization of the SK model decays as Universality 
strongly suggests that the same scaling should hold for 
our model. Although it seems not to be the case, at the 
first glance. Fig. [TJ] — top, a closer inspection confirms 
our expectation. Indeed, when plotted as a function of 
t~^/^ , see inset in Fig. [T2] — top, the thermoremanent 
magnetization curve has a finite non-vanishing slope at 
the origin. As we show in bottom panel of Fig. 1131 
finite size effects do not contradict this claim. In sum- 
mary, the magnetization decay for the hypercube suffers 
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Figure 12: Thermoremanent magnetization over Tlogf. The 
JANUS data (in red circles), follow a power law with an ex- 
ponent oc 1/r. Our results for D — 22 are shown in dark 
tonalities (lighter colors: _D<22). 




Figure 13: Thermoremanent magnetization at Tc vs. t~^^'^, 
for (up) different and D — 22, and (down) different sys- 
tem sizes for t„ = 4. The two insets are close-ups of the 
origin. 



from quite strong finite time effects, but asymptotically 
it scales with the proper exponent, at least at Tc. 



VII. NONEQUILIBRIUM CORRELATION 
FUNCTIONS AND FINITE SIZE EFFECTS 



The importance of finite size effects in nonequilibrium 
dynamics has been emphasized recently [23, HH . In our 
case, we have encountered important size effects, both in 
C{t,U), Fig. m and in ^{t^), Fig. [I3-top. 

We compare in Fig. [15] the finite D effects in C{t, iw) 
for two different MF models with fixed connectivity: the 
hypercube and a previously studied model (the random 
graph with connectivity z — 6, where each spin can inter- 
act with any other spin with uniform probability [29|). 
Clearly enough, the effects are much weaker in the hy- 
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Figure 14: Finite size effects in C{t,t„) at T = 0.7rc. 
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Figure 15: C{t,U) at T = 0.7Tc for t = U = 2^'^ vs. 1/L>^ 
We compare results obtained with two z = 6 models: one with 
hypercubic topology (red open circles) and another in a 
totally random graph (green full circles). 



percube model. 

It is interesting to point out that, although the finite 
size effects seems to be important in C{t,t^), they are 
largely absorbed when one eliminates the variable t in 
favor of C^{t,tvj), see Fig. [1^1 Hence, one of our main 
findings (the linear behavior of Cunk as function of C^) 
seems not endangered by finite size effects. 

A very clear finite size effect is in the coherence length, 
^(tw)- By definition, it cannot grow beyond D. Fur- 
thermore, what we find is that it hardly grows beyond 
D/2, Fig. [T71-top. Nevertheless, at short times, we 
can identify a _D-independent region, where it grows 
roughly as logtw Hence, one is tempted to conclude 
that ^£)=oo(^w) cx log^w At this point, finite size scaling 
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Figure 16: Ciink over C^(i,t„) at T = 0.7Tc for t„ = 2^^ and 
for different system sizes. 
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Figure 17: (Top) Coherence length ^ in the SG phase at 
T = 0.7rc vs. iw for different system sizes. (Bottom) same 
data of the top panel rescaled by D as a function of log t^/ D. 
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Figure 19: Fourier transform G{k) of C4{r) in units of the 
propagator [(p^ + 1/C^(^w)) ] [3 ^ £^ function of p, where 
the dimensionful momentum is p = k/D and ^(iw) ~ logtw. 
Recall that we are using postman metrics, hence, trans- 
lates to in the Euclidean metrics. We show results for two 
waiting times: t„ = 2^® (top) and t„ — 2'^'' (bottom). 



p of order I (k ^ D), although it seems to improve for 
smaller p as D grows. As far as we know, this is the first 
observation of the propagator in a numerical work. 



Figure 18: Fourier transform G(k) of C4{r) in units of G(0) 
as a function of \ogt^/D for several values of D and k at 
T = 0.7rc. For each value of fe, a different scaling function is 
found. 



suggests that both £,d/D and logty^/D arc dimcnsionlcss 
scaling variables. This is confirmed in Fig. [T71~bottom. 
where a spectacular data collapse occurs. This is fur- 
ther confirmed by the Fourier transform of Ci{r), G{k). 
Note that, since (74 (r) depends only on the length of the 
displacement vector r, also G{k) is rotationally invari- 
ant. Now, since k can range from to D, it is clearly a 
dimensionless quantity (a dimesionful momentum would 
be p = k/D). It follows that G{k)/G{0) is a dimension- 
less quantity that may depend only on a dimensionless 
variable, such as \ogt^/D. Our data support this expec- 
tation, see Fig. [THl 

As for the k dependence of G(fc), we expect a 1/p'^ 
behavior in the range of l/^(tw) P 1 (note 
that we are in the q = sector). Now, it is very impor- 
tant to recall that p'^ in Euclidean metrics translates into 
p^ in the postman metrics. We have also seen that the 
dimensionful p (postman metrics) corresponds to k/D. 
Thus, since in our range of tw, C(^w) ~ logiw, the prod- 
uct G{k) (jp' + 1/ log'^ tw) should be roughly constant as 
D grows. As we show in Fig. [121 the scaling is better for 



VIII. CONCLUSIONS 

We have studied a spin glass model in the D- 
dimensional unit hypercube in the limit of large _D, but 
with finite coordination number. We have shown that 
any short range model in such a lattice will behave as a 
mean field model in the thermodynamic limit (that coin- 
cides with the large D limit). An important advantage 
of this model is that it has a natural notion of spatial 
distance. 

We have shown that any statistical mechanics model 
on the hypercube with random connectivity would be 
afflicted by huge finite size effects, for purely geometrical 
reasons. The obvious cure has consisted in restricting 
the connectivity graphs to those with a fixed number of 
neighbors. Unfortunately, constructing such graphs is 
far from trivial. We have generated a subset of them by 
means of a simple dynamic Monte Carlo. In this way, we 
obtain graphs that are isotropic. We have checked that 
the Edwards- Anderson model defined over these finite 
connectivity hypercubes verify some consistency checks, 
including comparison with the analytically computable 
correlation function in the paramagnetic phase. 

We have numerically studied the nonequilibrium dy- 
namics in the spin glass phase. The three main features 
found were: (i) aging dynamics consists in the growth of 
a coherence length, much as in 3D systems, (ii) the seal- 
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ing of the two times correlation function implies infinitely 
many time-sectors, and (iii) the p** propagator has been 
observed. In addition, we have studied the finite size ef- 
fects in our model, finding that a naive finite size scaling 
ansatz accounts for our data. 



where K = tanh /3. In Bethe lattices, due to their cycle- 
less nature, Z'^^/Z'^ = 1 in the thermodynamic limit. 
Hence, we are left with the problem of counting the aver- 
age number of paths of length that join x and y, p{l.y). 
From it, we obtain 
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Appendix A: HIGH TEMPERATURE 
EXPANSION 

For sake of clarity, we will firstly discuss the calcu- 
lations for the random connectivity hypercube. Results 
for the fixed connectivity model will be then obtained by 
minor changes. 

Using the identity (/3 = l/T) 

^pj^ya^'yy ^ J^ya^ay tauh (3) , (Al) 

we can write the partition function and the spin propa- 
gator as: 



(cosh/?)^^ 



H n (1 + J^^'^^'^^ tanh/3) , (A2) 



{cr} (zw) 



CTxO'y = 



(A3) 

E{<t} Il{z^o) (1 + JzwCTzCTv, tanh /3) 



The high-temperature expansion (see, for instance [30|), 
expresses the propagator as a sum over lattice paths that 
join the points x and y, jx^y- 



(cTxCry) = Z-^ Z^J{ta.nh/3y 



(A4) 



where 1-^ represents the length of the path jx^y, J is the 
product of the coupHngs, Jzw, along the path, and is 
a restricted partition function obtained by summing only 
over all closed paths that do not have any common link 
with the path ^^^y 

However, when averaging over disorder, due to the ran- 
domness in the coupling signs, {(Tx(Jy) — 0. The spin 
glass propagator is obtained instead by averaging over 
disorder i^OxOy) . Clearly, the sum will be dominated 
by those diagrams where the go and return path are the 
same (thus, Jj^ = 1): 



{o^oyf = Z-^ ^7 hnh^/^]'^ = ^7' 

(A5) 



7a, -a 



C4(r) 



D 



(A6) 



l^>r 



The sum is restricted to l.y > r because the length of 
the shortest path that joins x and y is given by their 
postman distance r. 

In order to count the average number of paths, p{l^), 
let us distinguish two cases: — r and Ij > r. The first 
will give the leading contribution in the large D limit. 

The number of paths joining x and y in precisely r 
steps is r!, because the r steps are all taken along different 
directions and in a random order. For a given path, the 
probability of all the r links be active is {z/Dy . Hence 



(A7) 



Note that the D ^ factor compensates exactly the diver- 
gence of the (^) in Eq. (KE\i . 

In the case of 1^ > r, one has Ij = r + 2k, with k > 
0. Note that when — r the path contains r different 
directions (namely, the Euclidean components in which x 
and y differ). Each of these directions appear only once. 
However, when Ij > r, other directions must be included, 
we call them unnecessary. Note that, if the path is to 
end at the desired point, any unnecessary step must be 
undone later on. Hence, ^..^ — r is always an even number 
2k. Clearly, the number of such paths is bounded by 
r (r, k) , where F (r, k) is a ZJ-independent amplitude. 
On the other hand, the probability of finding all the links 
active is (zj DY^"^^ ■ Thus, we conclude that 



2k) = O 



( 1 



(AS) 



that results in a O (£> contribution to C/^(r). 
Then, in the large D limit we obtain (^A = zK): 



(74 (r) = A'- = e 



r loe A 



(A9) 



with finite size corrections of 0(Z?^^). Thus, we en- 
counter an exponential decay with an exponential corre- 
lation length given by 



cxp 



1 



|logA| 



(AlO) 



Summing all up, we can compute the spin-glass sus- 
ceptibility for the large D limit: 



oc oo 



r=0 



r=0 



-A' 



(All) 
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We see that when A = 1 the correlation no longer decays 
with distance, and the susceptibility diverges. Of course, 
one gets A = 1 precisely at the critical temperature, Tc, 
reported in Eq. 

The computation for the fixed connectivity model is 
very similar. One only needs to notice that, whereas 
the probability for the first link in a lattice path to be 
active is z/D, the probability for the next link is roughly 
{z — 1)/D (this is only accurate for large D). It follows 
that, again, the l-y — r paths are the only relevant paths 
in the high temperature expansion. We find that 



p{l^ = r) 



if r = Q, 
i^Y^^rl if r > 0. 



(A12) 



Again, we can use it to compute C4(r). In the large D 
limit, up to corrections of O (D^^), it is given by: 



1 if r = 0, 

^[{z-l)KY if r>0. 



(A13) 



which, taking A = {z — ^)K, also shows an exponential 
decay with 
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Using this spatial correlation function, we can either 
compute the SG-susceptibility in the fixed connectivity 
hypercube, 



r=0 



A 



z-l\-A 

or the integral correlation length, defined as ((T^ 



(A15) 



Er=oC4(r) 



X l~A 



(A16) 



Again, when A — 1. we find a critical point. The 
corresponding matches Eq. ([2]). The critical expo- 
nents, 7 = 1, V — 1, can be read directly from Eq. (|A15|) 
and (jA16[) . The reader might be puzzled by a mean field 
model with v ^ 1/2. The solution to the paradox is in 
our chosen metrics. Recall that the postman distance in 
the hypercube is the square of the Euclidean one. Hence, 
the correlation length in Eq. (|A16|) is the square of the 
Euclidean correlation length. 



Appendix B: SCALING AND DYNAMIC 
ULTRAMETRICITY 



As in Eq. let us assume that the spin time cor- 

relation function behaves for large as 



Figure 20: Parametric plot [x{t2),y{t2)] = 

[C(U-t2,t2),C{t2-t-i,t-i)], U > t2 > U with U fixed 
by the condition C(ti — tajts) = q and different tz- In the 
presence of dynamic ultrametricity, (|B2|l , the parametric plot 
should tend for large to the union oi x = q and y — q. 
The panels correspond to q = 0.25 (top, nice BB scaling 
but not ultrametric), q = 0.35 (middle, nice BB scaling and 
ultrametric) and q = 0.5 (bottom, supposedly ultrametric 
but poor BB scaling). 
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Figure 21: Dots: For each q and ta, as in Fig. 1201 we 
take the intercept with x = y, i.e. C* = C{ti — t2,t2) = 
C{t2 —t^jt^), and represent C* — 9 as a function of l/logta. 
Lines: analogous plot for the toy model described in the text, 
where the BB scaling is exact. 



C(i,tw) = /(a(t,iw)) , a(i,tw) =logt/logtw, (Bl) 
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Figure 22: For the data in Fig. 1211 we represent 
^(ti.t2)+o^(t2,i3) _ a{ti,t3) VS. l/logia. The dashed line cor- 
responds to Eq. (|B5|I . 



where the scaling function / is smooth and monotonically 
decreasing. From now on, we shall refer to this scaling 
as BB scaling (after Bertin-Bouchaud) . 

Let us see under which conditions BB scaling implies 
the ultrametricity property 

Cih~t3,t3) = min {Cih - h,h), C{t2 - t3, ^3)} , (B2) 

where ii 3> ^2 S> ^3 and ^3 tends to infinity. 

The natural time dependency is a power law choice 



h = h 

t2 - ^3 



(B3) 
(B4) 



with /ii > /i2. In that case, the large limit for the 
argument of the scaling function are: a{ti — t^^^t^i) = fii, 
ct{t2 - tsjts) = 112 and a{ti - t2,t2) = ^1 if /12 < 1 and 
ct{ti—t2, t2) = ^1/^2 if M2 > 1- Then, the condition (|B2p 
is only satisfied in case /Lt2 < 1- If, as it is the case for the 
critical trap model [11], f{a > 1) = constant [s^l, the BB 
scaling would imply dynamic ultrametricity. This is not 
the case for a general scaling function / such as, for in- 
stance, the one we get in Fig. [TOl Nevertheless, although 
this analysis implies that the dynamic ultrametricity is 
only present in our model in some range of parameters, 
let us try a more straight approach. 

We consider a fixed value for the correlation function, 
q. On the view of the previous considerations and of 
Fig. IIOI we should expect ultrametricity only for q > 
/(a = 1) 0.35. Now, for each ^3, we find ti such 
that C(ti — t^jt^) — q. Then, we perform a parametric 



plot of C{ti - t2,t2) vs. C{t2 - h,t3), for h < t2 < h. 
Ultrametricity predicts that, in the large limit, the 
curves should tend to a half square (e.g. the intersection 
of the straight lines x = q and y ~ q) and, in particular, 
when C{ti-t2,t2) = C{t2-h,t:i) ^ C* , C* should tend 
to q. 

We present in Fig. [20] results for three different val- 
ues of q: 0.5 (ultrametric region, but in our range of <w 
data do not scale according BB), 0.35 (ultrametric region 
and good BB scaling) and 0.25 (non ultrametric region 
but BB scaling works nicely). At the qualitative level, 
the parametric curves seem to tend to a corner, but the 
convergence is slow. Furthermore, there are no clear dif- 
ferences between the curves with q > f{a — 1) and those 
with q < f{a = I). Hence, we may try a more quantita- 
tive analysis. 

We obtain numerically C* , the point where C* = 
C[ti — t2,t2) = C(t2 — t3, t■i)^ and study C* — q as func- 
tion of 1/ log in Fig. [2TJ This choice is due to the fact 
that in the ultrametric region BB scaling predicts 



a{ti - t2, t2) = a{ti - t3,t3) + 



1 



21ogi3 



(B5) 



Hence, we expect that C* — q will be of order 1/ logt3 if 
ultrametricity holds. The numerical data confirms this 
picture only partly. For q — 0.35 the results are as ex- 
pected, yet for q — 0.25 the difference is decreasing fast 
as ts grows and it is hard to tell whether the extrapola- 
tion will be zero or not. For q — 0.5 (where BB scaling is 
not working for our numerical data) the behavior is non 
monotonic. 

To rationalize our finding, we consider a simplified 
model, where the BB scaling is supposed to hold exactly. 
The master curve /(a) is taken from the numerical data 
for C{t,t3 = 2^6) for D = 22. This toy model allows 
us consider ridiculously large values of ^3. As we see in 
Fig. [5TJ the peculiarities of the master curve cause a non 
monotonic behavior in q for an ample range of ^3 . 

The lack of monotonicity in q makes also on inter- 
est to focus on a, rather than on the correlation func- 
tion. With this aim, we consider the time ^2 where 
C{ti — t2,t2) ~ C{t2 — tsjts) = C*, and compute 
i [a(ti - t2, t2) + a{t2 - h.tz)] - a{ti - h^ts). BB scal- 
ing and ultrametricity combined, see Eq. (|B5p . imply 
that this quantity should be of order 1 / log ts (in the non 
ultrametric region, it should be of order one). Our results 
in Fig. [21] basically agree with these expectations. 
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